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1 . 0 SUMMARY 


As  alternatives  to  lengthy  series  expensions  for  globally  valid 
approximations  of  the  earth's  gravity  field,  piecewise  approximation 
methods  are  evaluated.  A single  global  equation  is  retained  only  for 
the  point  mass  and  second  zonal  terms  of  the  geopotential;  all  finer 
structure  undulations  are  modeled  by  a global  family  of  locally  valid 
functions.  A degree  23  spherical  harmonic  series  for  the  geopotential 
is  replaced  by  finite-element  approximations  within  the  spherical  shell 
out  to  1.2  earth  radii.  This  example  application  demonstrates  conclu- 
sively the  feasibility  and  desirability  of  tne  finite-element  approach. 

An  order  of  magnitude  reduction  in  the  calculation  time  for  gravitational 
acceleration  is  realized  over  conventional  calculations  with  spherical 
harmonic  recursions.  Preliminary  investigations  indicate  missile 
trajectories  and  satellite  orbits  can  be  integrated  efficiently  and 
accurately  using  finite-element  gravity  representations. 
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2 . 0 PREFACE 


This  report  constitutes  the  final  report  on  Phase  I of  Contract 
DAAG-53-76-C-0067  performed  by  the  University  of  Virginia  for  the  U.  S. 
Army  Engineer  Topographic  Laboratories,  Fort  Belvoir,  Virginia,  under  the 
sponsorship  of  the  Defense  Mapping  Agency.  The  Phase  I effort  represents 
the  culmination  of  a research  and  development  effort  initiated  under 
Contract  N00178-C00565  for  the  U.  S.  Naval  Surface  Weapon's  Center. 

The  Phase  I effort  is  concerned  with  replacing  global  gravity 
representations  (such  as  spherical  harmonic  series)  by  an  equivalent 
global  family  of  locally  valid  gravity  functions.  In  Phase  II  (documented 
in  a separate  report) , similar  methods  are  applied  to  modeling  the  fine 
structure  gravitational  anomalies.  In  both  the  macro  (Phase  I)  and 
micro  (Phase  II)  finite-element  modeling  efforts,  the  emphasis  has  been 
upon  formulating  gravity  fields  which  are  computationally  efficient. 

In  both  cases  gravity  models  have  been  developed  by  which  acceleration 
can  be  calculated  an  order  of  magnitude  more  quickly  than  by  traditional 
approaches,  without  significant  loss  of  accuracy. 

The  author  is  pleased  to  acknowledge  the  capable  technical  guidance 
provided  by  L.  A.  Gambino  of  the  U.  S.  Army  Engineering  Topographic 
Laboratory  who  served  as  technical  monitor  of  this  work.  The  computer 
programming  and  analysis  support  of  my  colleagues,  Mr.  J.  T.  Saunders 
and  Dr.  S.  Ray,  is  gratefully  acknowledged. 
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3.0  INTRODUCTION 


Gravitational  potential  of  an  arbitrary  body  cannot  be  written  as 
a closed  expression.  A number  of  infinite  series  have  been  formulated1 
and  truncations  of  such  series  currently  serve  as  models  of  the  gravity 
field  for  most  computational  purposes.  These  series  are  identical  to, 
or  are  motivated  by,  classical  product  solutions  of  Laplace's  equation 
in  spherical  coordinates1.  One  popular  form  of  the  spherical  harmonic 
series  for  the  geopotential  at  an  arbitrary  point  (r  = radius,  $ = 
geocentric  latitude,  X = east  longitude)  is 


N n 

0-?  I l If 

n=0  m=0  ' 


Pm  (sin  $)  [Cm  cos  mX  + sm  sin  mX]  (1) 

n n n 


where  R = a reference  value  of  the  earth’s  radius  = 6378.160  km,  GM  = 
the  earth's  gravitational-mass  constant  = 398601.2  km3/sec2,  P^(X)  = 
associated  legendre  functions  (non-no rmalized) , and  c”,  S™  = gravity 
coefficients  determined  to  render  simulated  satellite  motion  in  best 
agreement  with  observations. 

2 3 

Considerable  research  ' has  been  directed  toward  developing  stable 
recursion  algorithms  to  compute  Equation  (1)  and  the  south,  east,  and 
radial  acceleration  components  as 


„ 1 3U 

GS  = ' 7 3? 

G = 1 — 

E r cos  <j>  3X 


(2) 


In  spite  of  the  success  achieved  in  developing  feasible  algorithms 

based  upon  Equations  (1)  and  (2)  and  analogous  series,  the  increasing 

burden  of  computing  acceleration  from  ever  more  lengthy  global  gravity 

models  consumes  an  ever  larger  fraction  of  the  central  processor  time  in 

trajectory /orbit  integrations.  This  fact,  and  other  considerations  have 
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motivated  research  ' into  other  possible  global  gravity  representations 


The  present  report  departs  from  current  practice  by  separating  from 
the  onset  the  two  important  questions:  1)  What  gravity  representation 

should  be  fit  to  observed  data  to  estimate  unknown  model  constants  from 
satellite  observations  and/or  surface  gravimetry  (and  thereby  establish 
a global  gravity  field  approximation) ? 2)  In  highly  repetitive  gravita- 
tional calculations  (e.g.,  computing  acceleration  for  use  in  integration 
of  satellite  orbits/missile  trajectories) , what  gravity  model  is  most 
efficient  computationally? 

In  gravitational  modeling  research  to  date,  both  questions  have 
been  (perhaps  subconsciously)  addressed  simultaneously.  Explicit  sepa- 
ration of  the  "optimum  determination"  and  the  "optimum  use"  questions 
appears  to  be  a most  important  consideration  in  further  refinement  of 
existing  gravity  models  to  accommodate  the  ever  more  precise  and  abundant 
observed  data.  Without  engaging  in  the  important  quest  for  a model  which 
can  be  best  determined  from  observed  data,  we  address  the  problem  of 
determining  an  "optimum  use"  gravitation  model. 


4.0  FINITE-ELEMENT  APPROACH 

The  dominant  global  features  of  the  gravity  field  are  efficiently 
represented  by  the  dominant  low-degree  terms  in  Equation  (1) . This  fact 
motivates  the  segregation  of  the  total  gravitational  potential  at  a point 
(r,  $,  X)  as 

U = Uggptr,  <t>,  X)  + AU(r,  <f>,  X)  (3) 

where  it  is  understood  that  a single  low-degree  truncation  of  Equation 

(1)  is  adopted  to  define  U__„  globally.  But  instead  of  attempting  to 

model  "everything  else"  as  a single  global  series,  we  shall  determine 

a global  family  of  locally  valid  disturbance  functions  to  model  AU  and 

its  gradient.  Since  the  local  equations  must  model  only  the  gravity 

undulations  (in  addition  to  U___)  in  a specific  local  volume,  it  is 

REF 

reasonable  to  anticipate  significantly  more  compact  expressions  than  a 
single  global  expansion  of  comparable  local  accuracy.  This  is,  of  course, 
the  thesis  of  the  finite-element  approach. 

Having  decided  to  pursue  the  finite-element  approach,  it  is  necessary 
to  make  several  important  decisions;  it  is  necessary  to  define  1)  what 
portion  of  the  geopotential  is  to  be  approximated  as  AU  in  Equation  (3) , 

2)  the  specific  mathematical  structure  to  use  as  local  approximations  of 
AU  and  its  gradient,  3)  the  size  of  the  finite  elements,  4)  the  procedure 
to  be  used  in  determining  numerical  values  for  coefficients  in  the  local 
approximations,  5)  the  order  of  continuity  requirements  between  adjacent 
approximations  across  their  mutual  boundaries  of  validity,  and  6) 
procedures  to  evaluate  the  validity  of  the  finite-element  model. 

These  decisions  are  coupled  and  affect  the  accuracy  and  efficiency 
of  the  resulting  finite-element  model  in  a complicated  way.  Analytical 
attempts  to  resolve  these  issues  a priori  proved  unsuccessful.  Therefore, 
guided  by  intuition,  we  have  developed  experimental  software  and,  based 
upon  systematically  collected  empirical  data,  we  arrived  at  a prototype 
finite-element  model  of  the  geopotential.  We  present  several  approaches 
for  constructing  gravitational  finite-element  models  and  summarize 
numerical  experiments  with  them  in  the  following. 
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5.0  WEIGHTING  FUNCTION  APPROACH  TO  PIECEWISE  CONTINUOUS  APPROXIMATION 


Recent  papers  document  the  theoretical  development  and  application 
of  a versatile  piecewise  continuous  approximation  technique.  This 
weighting  function  approach  determines  an  arbitrarily  large  family  of 
locally  valid  functions  which  join  with  rigorous  piecewise  continuity 
through  any  prescribed  order  of  partial  differentiation.  As  applied  to 
the  approximation  of  a three  variable  function,  FfX^  X2>  X^)  each  final 
local  approximation  F(X^,  X2,  X^)  is  defined  as  a weighted  average  of 


eight  preliminary  local  approximations  {f^  (x^, 


X2,  X3);  i,j,k,  = 0,1} 


FIX  X X , . I I J W ,X  ,X  ,X  | 

1=0  3=0  k=0  J 


£iJk<Xl'X2-X3>  141 


The  eight  preliminary  approximations  may  have  arbitrary  mathematical 

structure  (a  most  flexible  and  powerful  feature  of  this  approach)  but 

must  be  determined  in  such  a fashion  that  their  respective  centroids  of 

validity  lie  at  the  eight  vertices  of  a parallelopiped  (in  X^,X2,X3, space) , 

which  defines  the  volume  in  which  the  final  local  approximation  ?(X  , X ,X  ) 

6 7 x z j 

is  valid.  The  form  of  the  weight  functions  can  be  selected  ' to 
guarantee  that  F and  its  similarly  determined  adjacent  approximations 
are  continuous  across  their  mutual  boundaries  of  validity;  the  order  of 
partial  differentiation  to  which  continuity  is  desired  is  controlled  by 
selecting  appropriate  weight  functions.  Specifically,  Table  1 gives 
the  weight  function  W^^fX^X^X^)  for  various  orders  of  continuity.* 

The  remaining  seven  weight  functions  are  obtained  by  reflecting 

W111(X1'  X2'  X35  as 


WYW  = wui(1  - V 1 - X2 ' 1 " V 


WVW  = Wlll(1  - V 1 - X2 ' V 


WYW  - wm(1 " xi'  V 1 " V 


tUnder  the  assumption  that  X^,  Xj,  X^  are  nondimensional  local  coordinates 
with  (X^  = i,  X2  = j,  X = k)  locating  the  eight  vertices  of  the  unit 
cube  of  validity  of  F(X^,  X2,  X3) . 


“ioo'VW  ■ 1 - V 1 - V 


(Sd) 


'WVW  - "m(1  - *1'  V V 


"lOl'VW  ■ “lll(V  1 - X2'  X3> 


“U0«X1'X2'X3>  * "iu(Xl'  X2'  1 " V (59) 

The  weight  functions  are  positive  and  satisfy  the  constraint  that 
111 

I l l wiik<xi'  V X3}  = 1 (6) 

i=0  j=0  k=0  13 K 1 2 3 

They  may  be  interpreted  geometrically  as  follows:  The  maximum  value  W^_.^ 

is’  unity  which  occurs  at  X,  = i,  X2  = 3 , X3  = k;  the  surfaces  of  constant 
weight  are  spherical  in  the  vicinity  of  (i,j,k)  but  become  increasingly 
angular  until  the  surface  of  zero  weight  is  the  walls  of  the  cube 
opposite  to  i,j,k.  Thus,  in  Equation  (4),  W.  causes  the  preliminary 
local  approximation  f„^  to  dominate  F in  the  vicinity  of  f ^’s  centroid 
of  validity,  but  has  no  effect  on  F (in  value  or  first  m partial  deri- 
vatives) along  the  opposite  cell  boundary.  The  key  to  piecewise  continu- 
ity is  the  fact  that  F is  completely  defined  (on  each  of  the  six  boundary 
"walls")  by  the  four  preliminary  approximations  whose  centroids  are  the 
vertices  of  the  respective  walls.  Since  the  four  preliminary  approxi- 
mations (defining  any  given  cell  wall)  are  shared  by  adjacent  final 
approximations,  it  is  clear  that  piecewise  continuity  is  assured. 


6.0  LOCAL  GRAVITY  APPROXIMATIONS  VIA  TAYLOR'S  SERIES 


Given  an  a priori-determined  global  gravity  model  (e.g.,  a spherical 
harmonic  series) , perhaps  the  most  obvious  strategy  of  generating  local 
approximations  is  by  Taylor’s  series.  The  disturbance  potential  and  its 
gradient  can  be  locally  approximated  as  the  truncated  Taylor's  series 
of  three  variables 


4=ijk(r,*,X>  = ,y  ♦ l I l gijkxi  x2  x3  <7> 


n=l  1=1  J=1 


where 


AG. 


jk 


r au  i 

f Au  ^ 

-l/r  (3  (AU)/3  <j>) 

> = < 

AG. 

s 

1/r  cos  <t>  (3(AU)/3  X) 

r 

age  f 

3 (AU)/3r 

AG_ 

J 

RJ 

rAU  ^ 

(Ar)I(A4i)  J(AX)K  3n 

AGC 

(8) 


...JK 


I ! J! K! 


o'rI3(j)J3XK  ^ AGr 
AG, 


(9) 


R J (r , (}> , A ) = (ri , <(> j /Xk) 


K = n-  I - J,  M = order  of  the  local  Taylor’s  series.  (r.,$.,X ) = an 
arbitrary  local  expansion  point  [which  is  the  centroid  of  validity  of 
Equation  (7)],  (2Ar , 2A<£ , 2AX)  = dimensions  of  the  region  of  validity  of 
Equation  (7),  and  (X^X^X^)  = {(r  - r^/Ar,  (<j>  - <p.)/A<fi,  (X  - X^J/AX}  = 
nondimens ional  local  coordinates. 

The  partial  derivatives  Equation  (9)  can  be  regorously  derived  from 
the  parent  global  model  of  ACJ.  Reference  14  gives  analytical  expressions 
for  computation  of  the  partial  derivatives  specifically  for  the  case 
that  the  parent  global  model  of  AU  is  a spherical  harmonic  series.  It 
should  be  pointed  out  that  specific  numerical  values  for  the  elements 
of  G__.  are  computed  a priori  and  stored;  for  centroids  of  validity 
distributed  over  the  (r,4>,X)  space  according  to  some  specified  pattern. 

In  using  (say)  the  foregoing  weighting  function  formulation  to 
compute  local  disturbance  acceleration  (from  a piecewise  continuous, 
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finite-element  model) , the  appropriate,  previously  computed,  set  of 
coefficients  are  employed  to  compute  Equation  (7)  as  preliminary 
local  approximations  for  substitution  into 


r U(r,<(.,X)> 

u^lr.+  .A)  ' 

G (r  ,4>,  X) 

> = l 

-1/r(3tW3*>  | 

GE(r,<fr,X) 

f S 

1/r  cos  ((•OU^/aX)  | 

J3R(r,<fr,X)^ 

v.  REF  J 

l I I wiik 

i=0  j=0  k=0  J 


(X1,X2,X3)AGijk(r,4.,X)  (10) 

where 

= (r  - rQ)/Ar,  X2  = (<*>  - 4>Q) /A4> , X3  = (X  - XQ)/AX 

are  nondimensional  local  coordinates  and  ^q'^q'^q)  are  coordinates  of 
the  "lower  left  corner"  of  Equation  (l)’s  region  of  validity:  {0  < X.  < 1; 
i = 1,2,3}. 

The  gravity  representation  (10)  leads  to  a nonuniform  distribution 
of  errors  over  the  cell  volume.  Observe  that  the  approximations  become 
exact  as  the  displacement  of  the  evaluation  point  from  a centroid  of 
validity  (expansion  point)  decreases  to  zero;  but  more  generally,  the 
final  approximation  is  the  average  of  eight  approximations  containing 
errors.  This  observation  led  us  to  expect  from  the  outset  that  Taylor's 
series  would  not  be  the  optimum  choice  of  preliminary  approximation 
functions;  one  should  seek  preliminary  approximations  with  more  uniform 
error  distributions. 
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7.0  LOCAL  GRAVITY  APPROXIMATION  VIA  LEAST  SQUARES  APPROXIMATIONS 


As  an  alternative  to  the  local  Taylor's  series  [Equation  (7)],  we 
consider  as  the  local  model  of  disturbance  gravity 


fu  > 

r u 

IJK 

M n n-I 

- [ l I J 

SIJK 

n=0  1=0  J=0  ' 

eijk 

^gr> 

^RiJx> 

PIJK(WV 


(11) 


where  K = n - I - J. 

[F000 ; F001 ' F010 ' F1000 ; F002 ' F011  ,F101 ' F020 ' F110 ' F200 FM00 ( X! ' X2  ' X3 ] 1 
are  a suitable  set  of  linearly  independent  basis  functions,  and 

= tU000U001U010U100" 'UM00]  ' 

= [S000S001S010S100‘ ‘ *SM00] ' 

= CEoooEooiEoioEioo‘ ‘ ^MOO1 ' 

{r}  = [R„„„R„„ , R. . _R.  _ . . . . R,„ . - ] , 

000  001  010  100  M00 

are  coefficients  determined  so  that  the  sum  square  residual  error  (between 
Equation  (11)  and  the  parent  model  of  disturbance  gravity)  is  minimized. 

In  particular,  if  the  least  square  coefficient  estimates  are 

determined  by  fitting  Equation  (11)  to  local  evaluations  of  a global 

gravity  model,  then  the  coefficient  estimates  are  given  by  the  normal 
11 

equations 


{U}  = [B]{AUCJ 
{S}  = [B]{AGSC} 
{E}  = [B]{AGEC> 
{R}  = [B]{AGRC} 


(12) 


where 


[B] 


T -IT 
[(A  WA)  A W] 


(13) 
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F000(Xll'X21,X3r  F100(X11'X21'X31)  *“  FM00w‘11'“21'“31' 

F000(X12'X22'X32)  F100 (X12 'X22 ' V ***  FM00 <X12'X22'X32> 


[A]  = 


F000(Xln'X2n'X3n)  F100 (Xln'X2n'X3n)  FM00 (Xln'X2n'X3n) 


[W]  = [w  ] = an  n x n positive  definite  weight  matrix. 

(S0C>T  . WUUU,X21.X31)  4U(X12,X22,X32)...  4U(Xln,X2n,X2n))' 

{4GSc)T  - UOs(Xu,X21X  31)  4GS(X12,X22,X32)...  AGs(Xln,X2n,X3n)} 

T ) 

W<3V  ‘ (4VXU'X21'X31>  4VX12'X22'X32>-"  4VXl,,'X2n'X3nl  1 

{AGE  }T  - (AGR(X  ,X  ,X3  ) AG  (X12.X22,X  ) A=R(Xl„’X2n'X3n> 1 

J 

are  local  evaluations  of  disturbance  potential  and  acceleration  (from  a 
parent,  global  gravity  model)  at  the  set  of  points 


<XL£'X2-£'X3.£>  = 


:l  ~ rj  *1  ~ Xl  ~ Xk 
dr  d<$  dA 


; •£  = 1,2, ...  ,n  (16) 


in  the  vicinity  of  the  centroid  of  validity  (r^,$.,A^).  It  has  been 
found  advantageous  to  generate  the  data  (15)  on  a uniform  grid  in  the 
local  (X^£,X 2£'X3£)  space;  providing  this  grid  is  held  identical  for  all 
cells,  then  Equations  (13)  and  (14)  can  be  computed  only  once  and  simply 
reused  in  (12)  to  operate  upon  appropriate  local  data  to  generate  the 
entire  global  set  of  local  coefficients  for  use  in  Equation  (11) . 

An  infinity  of  choices  exist  for  the  basis  functions  in  (11) ; after 
some  experimentation,  we  selected  the  set  of  all  Chebyshev  polynomial 
products  up  to  Mth  order  (Table  2)  as 


VX3> 

vv 

w 


j 


: 


i ■- 

i 

' ' 


' i 

I I 


t (X) 
n 


(dt  /dX) 
n 


2X2  - 1 
4X3  - 3X 


4X 


12X2  - 3 


Recursions : 


t (X)  = 2X  t , (X)  - t _ (X) 
n n-l  n-<£ 


•1  < x < 1,  n > 2 


dt  dt  . 

(dt  /dX)  = 2t  (X)  + 2X  — ~ n > 3 

n n-l  ax  dX 


Shifted  Chebyshev  polynomials: 


T (Y)  = t (2Y  - 1)  0 < Y < 1 

n n - - 


Table  2 Chebyshev  Polynomials. 
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8.0  NUMERICAL  TRADEOFF  STUDIES 


To  gather  the  empirical  evidence  upon  which  to  base  selections  from 
the  many  alternatives  implicit  in  the  above  developments,  a versatile 
computer  software  systemt  has  been  developed.  Figure  1 summarizes  a 
portion  of  a tradeoff  study  whose  objective  is  to  help  decide  whether 
Taylor's  series  or  locally  fit  Chebyshev  polynomials  should  be  employed 
as  local  gravity  approximations.  The  RMS  of  the  acceleration  error  norm 
was  determined  by  computing  the  variance  of  a directly  evaluated  accelera- 
tion error  sample.  In  this  case,  the  sample  consisted  of  216  uniformly 
spaced  error  calculation  points  over  each  finite  element.  Using  the 
RMS  acceleration  error  norm  as  an  accuracy  criterion  (for  a given  order 
of  the  local  approximation)  it  is  clear  that  the  locally  fit  Chebyshev 
polynomials  are  superior  to  equal  order  Taylor's  series.  The  very 
nonuniform  error  distribution  characteristic  of  truncated  Taylor's  series 
(zero  error  at  expansion  point,  but  rapidly  degrading  away  from  that 
point)  is  partially  compensated  for  in  the  weighting  function  approach 
by  the  bell  shaped  weights  and  redundancy  of  the  method.  However,  the 
RMS  of  the  locally  fit  Chebyshev  polynomials  have  been  found  consistently 
superior  for  low  order  approximations  (<  4) , usually  by  an  order  of 
magnitude.  Consideration  of  other  goodness  of  fit  criteria  (smallness  of 
mean  error,  smallness  of  maximum  error,  near  Gaussian  residuals,  etc.), 
support  the  choice  of  locally  fit  polynomials  over  local  Taylor's  series. 

Based  upon  the  numerical  experiments  done  thus  far,  the  following 
tentative  conclusions  have  been  drawn: 

1)  Rigorous  piecewise  continuity  of  the  local  approximations, 
while  desireable  from  conceptual  and  esthetic  viewpoints,  appears 
weakly  justified  using  the  small  RMS  error  criterion.  The  weighted 
average  of  eight  local  approximations  is  consistently  superior  to 
any  of  the  original  approximations,  but  often  not  sufficiently 
superior  to  justify  the  8:1  redundancy  of  the  method. 


tThe  FORTRAN  IV  language  has  been  used  exclusively,  all  computation  has 
been  carried  out  using  the  University  of  Virginia's  Control  Data 
Corporation  CYBER  172  computer  system. 
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[Cell  Size:  5'°x5°x0.2R  Located  on  the  Equator] 


Figure  1 Order  of  Local  Approximations  vs  RMS  Acceleration 

Approximation  Error,  e = (spherical  harmonic  model) 
minus  (finite-element  acceleration  model) . 


V 


T 


2)  By  fixing  the  order  of  the  local  approximations  and  selecting 
a specific  accuracy  criterion,  straightforward  variation  in  the 
finite-element  dimensions  leads  quickly  to  the  family  of  maximum 
volume  elements  consistent  with  these  two  constraints.  Global 
numerical  tests  support  the  conclusion  that  global  decisions  can 
usually  be  reliably  made  based  upon  local  numerical  experiments,  so 
long  as  the  local  experiments  include  a full  range  of  latitude 
variation. 

3)  Analysis  of  repeated  statistical  samples  of  acceleration  errors, 

(between  local  Chebyshev  approximations  and  a degree  23  global 
spherical  harmonic  series)  taken  from  various  finite  elements  support 
the  following  conclusions  regarding  distribution  of  approximation 
errors:  (a)  The  mean  acceleration  error  has  always  been  found  to 

be  at  least  one  order  of  magnitude  less  than  the  sample  RMS  accelera- 
tion error,  providing  the  sample  was  taken  from  at  least  200  sample 
points  located  either  on  a uniform  grid  or  randomly  located  within 
the  element;  and  (b)  The  maximum  acceleration  error  has  always  been 
found  to  be  less  than  4 times  the  sample  RMS  acceleration  error, 
with  the  same  sample  restrictions  as  (3a). 
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9.0  PROTOTYPE  FINITE-ELEMENT  MODEL  OF  THE  GEOPOTENTIAL 

Based  upon  insight  gained  in  parametric  studies  with  the  experimental 
software,  a prototype  finite-element  model  has  been  developed.  The  global 
finite  model  has  the  following  structure: 


r U(r,4>,X)  "N 

r U(r,*,Xp 

ru 

IJK 

Gs(r ,*,X) 

r< 

Gs  (r,4>  ,X) 

M n n-1 

U i i 1 1 

SIJK 

Gg(r,*,X) 

Ge  ( r , , X ) 

' n=0  1=0  J=0 

eijk 

G (r,<j>,X) 

V.  R J 

G (r,*,X) 

V R J 

REF 

^ Rjjk> 

WWW 


(18) 


where 


uBEP(r**»*>  s T [1  + C2  7 P“(sin4,)] 


SREF 

^REF 

GRrEF 


1 3UREF  _ H3M  0 

r 3*  ' r2  2 


7/ 


3U 


REF 


r cos  <(>  3A 
3U. 


= 0 


REF  GM  „ ,„0 

t—  1 + 3C. 

2 2 


(sin  <(»)  ] 


(19) 


(20) 


(21) 


(22) 


and  U , S , E , R^  are  the  appropriate  set  of  a priori  computed 
X J K IJK  ijk  xjjn. 

local  coefficients;  computed  according  to  Equations  (12-14)  with 


WVW  = WWW 


,12 


to  accurately  replace  the  spherical  harmonic  model  of  disturbance 
gravity 


23  n 


l l 

n=2  m=l 


_ » n 

— P (sin  $) [C  cos  mA  + S sin  mX] 
r ) n n n 

(24) 


and  its  gradient. 


Figure  2 is  a projection  of  a global  contour  map  of  the  radial 
disturbance  acceleration  [3(AU)/3r]  on  the  earth's  surface.  To  fully 
define  the  procedure  for  constructing  the  finite  elements,  the  following 
decisions  were  made:  1)  develop  a finite-element  model  for  the  spherical 
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shell  within  1.2  earth  radii;  2)  adjust  the  finite-element  size  and/or 
the  order  of  the  local  approximation  so  that  the  acceleration  approxi- 
mation errors  enter  (at  worst)  in  the  seventh  significant  figure;  and 
3)  fix  the  order  of  the  local  approximations  at  M = 3 [in  Equation  (19)]. 

Holding  the  radial  dimension  fixed  at  0.2  R and  adjusting  the 
longitude  by  latitude  dimensions  to  maintain  requirement  (2)  led  to  the 
set  of  1500  finite-elements  whose  bases  are  shown  in  the  flat  projection 
of  Figure  3.  The  acceleration  error  residual  norm 

AG  = [ (GSH  - GFE)T(GSH  - GFE) ]%  (25) 


where 


GSH  = 


spherical  harmonic  series  model 


GFE  = 


finite-element  model 


were  computed  as  were  its  RMS  and  maximum  value  for  N = 23  in  Equations 
(1)  and  (2) . We  found  that 


. n 

AG  = i 7 AG.  = 0.00000003  m/sec2 

N . X 


RMS  = 
AG 


" N 

= i 7 AG.2 
N i 

x=l 


= 0.000002  m/sec2 


AG„-„  = 0.000008  m/sec2 

MAX 

N = 360  x 91  = 32760  sample  points  (l°xl°  grid  in  the 

Northern  Hemisphere) 

These  errors  are  the  worst  case  errors  arising  generally  in  the  most 
anomalous  region  near  the  earth's  surface;  analogous  statistical  analyses 
reveal  the  magnitude  of  the  gravity  modeling  errors  on  the  surface  of  a 
1.03  R sphere  are  about  one  half  of  the  above  values  while  the  errors 
are  reduced  by  order  of  magnitude  on  the  surface  of  the  1.2  R sphere. 
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10 . 0 TRAJECTORY/ORBIT  INTEGRATION 


This  level  of  precision  has  been  found  to  be  satisfactory  for  inte- 
gration of  most  missile  trajectories  and  satellite  orbits;  typical  single 
revolution  position  integration  errors  have  been  found  to  be  on  the  order 
of  0-25  m.  Clearly,  greater  precision  can  be  easily  achieved,  if  desired, 
by  decreasing  the  finite-element  size  (thereby  increasing  the  number  of 
elements) , or  by  increasing  the  order  of  local  approximations  (thereby 
decreasing  the  computational  efficiency  of  the  method) . 

The  computational  speed  of  the  finite-element  model  versus  typical 
spherical  harmonic  recursion3  favors  the  finite-element  model  of  Equation 
(19)  by  better  than  an  order  of  magnitude.  The  computational  picture 
is  complicated,  however,  by  virture  of  the  fact  that  random  access 
retrieval  of  previously  stored  coefficient  subsets  is  necessary.  In 
this  case,  each  component  of  acceleration  requires  a total  of  (20) (1500)  = 
30,000  coefficients  to  define  the  entire  global  family  of  gravity  functions 
(although  only  twenty  are  used  in  each  element) . Since  the  elements  are 
large  (hundreds  of  miles)  compared  to  small  errors  (tens  of  miles) 
associated  with  two  body  or  other  simplified  dynamic  extrapolations, 
simple  logic  can  be  devised  to  bring  into  central  memory  several  local 
sets  of  coefficients  before  they  are  needed  and  thereby  hold  the  lost 
time  during  random  access  to  a minimum.  For  the  ballistic  missile  problem, 
most  of  the  acceleration  evaluations  occur  in  two  local  regions  (i.e., 
during  atmospheric  powered  flight,  and  during  re-entry) ; thus  hundreds 
of  acceleration  evaluations  are  likely  within  a single  finite  element. 

Even  if  no  systematic  pre-access  scheme  is  used  and  the  local  coefficients 
are  accessed  each  time  acceleration  is  required,  the  finite-element 
approach  still  maintains  an  order  of  magnitude  advantage  on  the  UVA  CYBER 
172  computer  system. 
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11.0  CONCLUSIONS 


The  finite-element  approach  to  modeling  the  geopotential  has  been 
studied  analytically  and  numerically.  Many  degrees  of  freedom  exist  in 
our  approach  to  this  problem.  The  specific  finite-element  model  developed 
and  discussed  herein  is  not  put  forth  as  the  optimal  computational  model 
of  the  geopotential.  Rather,  we  believe  the  prototype  finite  model  to 
be  a representative  finite-element  geopotential  model  which  will  probably 
be  improved  upon  in  future  refinements  of  our  approach.  The  computational 
speed  advantages  over  spherical  harmonic  expansions  is  clear,  however 
(primarily  because  a 23rd  order  expansion  is  locally  replaced  by  a 3rd 
order  expansion) . The  ultimate  computational  speed  advantage  depends 
upon  the  number  of  random  accesses  required  to  maintain  the  appropriate 
local  coefficients  in  core  memory.  For  a given  trajectory/orbit  integra- 
tion, simulations  done  to  date  support  the  conclusion  that  an  order  of 
magnitude  savings  is  achieved. 
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